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About 20% of all massive stars in the Milky Way have unusually high ve- 
locities, the origin of which has puzzled astronomers for half a century. 
We argue that these velocities originate from strong gravitational interac- 
tions between single stars and binaries in the centers of star clusters. The 
ejecting binary forms naturally during the collapse of a young ( ^ 1 Myr) 
star cluster. This model replicates the key characteristics of OB runaways 
in our galaxy and it explains the 100 runaway stars around young 
star clusters, e.g. R136 and Westerlund 2. The high proportion and the 
distributions in mass and velocity of runaways in the Milky Way is re- 
produced if the majority of massive stars are born in dense and relatively 
low-mass (5000 - 10000 M©) clusters. 

Most stars in our galaxy have a relatively low velocity. However, there is a population of 
fast moving stars, called OB runaways (7), which have considerably higher space motions 
of > 30km/s (2). The origin of such velocities can be attained in two very distinct ways: 
A runaway can be launched when its binary companion explodes in a supernova (3), or by 
the ejection via a dynamical slingshot (4). The relative importance of both mechanisms has 
remained elusive, mainly because both are associated with young stellar populations and 
the high speed of a star is generally observed long after it has moved away from its birth 
place. 

A massive star can be accelerated effectively by a three-body dynamical interaction (4). 
Every hard interaction eventually results in a collision between two or all three participating 
stars, or in the escape of one star and one binary (5). The velocity acquired by the ejected 
star easily exceeds the escape speed of a star cluster: For a binary with total mass M5 and 
semi-major axis a, the typical velocity with which the single star is ejected is v'^- = GM5 / a. 

It is typically the least massive star that is ejected (5), and both components of the 
ejecting (bully) binary are therefore likely to be more massive than the escaping star. For 
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higher-mass binaries, the maximum cross section for ejecting an intruder with high speed 
shifts to a larger orbital period, and to a higher binding energy (Supplement § A). The 
binary that most effectively produces massive runaways is relatively wide, composed of 
the most massive stars in the cluster and located in the core. The gravothermal collapse of 
a cluster core naturally produces such a binary (6), which subsequently hardens by ejecting 
stars (S), until it experiences a collision or hardens to ^ 10, 000 kT (7), upon which an 
encounter causes its ejection from the cluster (9) 

A binary with a binding energy of 10 to 10,000 kT is therefore a natural consequence 
of the core collapse of a star cluster (70), and leads to the frequent ejection of other cluster 
members, until it eventually ejects itself or engages in a phase of Roche-lobe overflow. 
During the hardening process, each encounter typically drains ^ 40% (5) of the binary's 
binding energy. In this process it ejects on the order of log(10, OOOkT/lOkT)/ log(l + 
0.4) 21 stars from the cluster before the binary ejects itself or becomes semi-detached. 
The collapsed cluster core can support only one bully binary at a time and as a consequence 
the number of runaways produced does not depend on cluster mass. 

We quantized these results by performing a series of A^-body simulations for clusters 
of 6,300 Mq to 2.0 X 10^ Mq (Tab. SI). For each simulation, we recorded the moment 
of core collapse and the time of ejection, mass and velocity of the escaping stars. Each 
cluster produced a single binary that ejected ^ 23.0 ± 5.9 stars with ^ej > 30km/s of 
which ^ 5.2 ± 1.6 stars have a mass m > SMq, irrespective of the mass of the cluster 
(Supplement § B); this number is consistent with the 6 runaways with m > 8 Mq ejected 
from R136 (11, 12). With a mass of M 60, 000 Mq and a core radius of rcore ^ 0.10 pc 
{13, 14) R136 represents the population of young (~ 2 Myr), dense and massive star clusters 
in the local group of galaxies (Fig. S2). This star cluster is ideally suited for comparing with 
our simulations because of the excellent quality of the copious observations. 

In our simulations, the large mass and relatively wide orbits of the binaries produced 
during the gravothermal collapse of a cluster core completely dominate the cross section 
for stellar encounters (Supplement § A). The encounter probability then should depend only 
weakly on the mass and size of the single stars in the cluster core. As a result, the mass 
distribution of the ejected stars roughly follow the mass function in the cluster core which 
should be relatively flat compared to the global cluster mass function (75), because equipar- 
tition of kinetic energy has resulted in the massive stars to sink to the cluster center (Sup- 
plement §B). 

In Fig. [H we compare the distribution of the masses of the runaway stars with the mass 
function in the cluster core (see also Supplement § B). The mass distribution of the runaway 
stars is statistically indistinguishable from the core distribution. For stars of m 10 the 
mass distribution of the runaways are flatter than the Salpeter slope, whereas it is consid- 
erably steeper for m ^ 30 Mq. The steep high-mass end is consistent with the recently 
observed distribution of massive runaway stars and field stars in the small Magellanic 
cloud (16). 

The binaries from the A^-body simulations populate the area around the dashed curve 
in Fig. El These binaries came from a longer orbital period and lower mass to reach their 
minimum orbital period mediated by encounters and collisions (Supplement § A). Binaries 
to the right of this curve are in the process of dynamical hardening by ejecting stars from 
the cluster; they most efficiently produce runaways. When they move away from the dashed 
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Figure 1: The cumulative mass functions (1 — F{m)) of runaway and core (r < O.lpc) 
stellar populations. The solid black lines indicates the Salpeter initial mass function (left) 
and the observed mass function for massive field stars and runaways in the small Magellanic 
cloud (7(5) (right black curve). The high-mass end of the function of the runaway stars 
from the simulations are consistent with the core mass function (dotted red curve); all are 
strongly mass segregated. The break point of the power law for the runaway mass functions 
is around IOMq and agrees with the expectation from the energy equipartition (77, 18). 
The green dash-dotted, red solid and blue dashed curves give the results for simulated star 
cluster with a total mass of 6 300 (model A), 5. 1 x 10^ (model C) and 2.0 x 10^ 
(model D) respectively (Tab. SI). 
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Figure 2: The mass and orbital period of the binaries observed in nearby massive star 
clusters. The observed binaries are from the star clusters R136 (red squares) (79-22), 
Trumpler 16 (blue up-pointed triangles) {23), NGC3603 (magenta bullet) (20), Wester- 
lund 1 (yellow down-pointed triangles) (24), Westerlund 2 (green up-pointed triangle) (25) 
and Trumpler 14 (green star) (26). (For the binary RMS we have plotted also the latest 
preliminary parameters (27) as the red right-pointed triangle.) The simulation results are 
represented with the letters A, B, C and D for the respective models (Tab. SI). The arrow 
(indicated with M) gives the change in mass and orbital period of a binary as a result of 
adiabatic mass loss by a stellar wind. The blue diagonal dashed line (right) indicates the 
parameters for which the cross section for ejecting a single star with ^gj > 30km/s is maxi- 
mal, adopting an equal-mass interaction (Supplement § A). The solid black curves give the 
encounter rate (as indicated) for parameters in the core of R136 (adopting a core radius of 
0.1 pc). 
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blue curve their encounter frequency drops and stellar mass loss in the winds of the stars 
start to dominate their evolution towards lower masses and longer orbital periods (Fig.O. 

The short ^ 10 day period binaries observed inside young star clusters like R136, 
must be primordial and probably never experienced a dynamical encounter with another 
cluster member (Fig. [2]). From a dynamical perspective these binaries can be considered as 
single stars (70), because their interaction cross section is very small (Fig. SI). A strong 
interaction with such a tight binary tends to cause a collision between two or all three 
stars (28, 29); ^ 80% of the collision products remain in the cluster, the others escape 
(Supplement § A). 

The relatively wide (P ^ 90 year) binary HD93129AB (26, 30) (green star in Fig.© 
probably formed during the early collapse of Trumpler 14 (31, 32), With the observed mass 
and density of Trumpler 14 (33) this binary would be about 2,700 kT, consistent with being 
formed dynamically during core collapse and currently in the process of hardening, but 
insufficiently hard to eject itself from the cluster. We expect that HD93129AB has ejected 
^ 16 stars, a few of which with m > 8 Mq. 

The outskirts of the star cluster R136 also exhibits evidence of its violent dynami- 
cal history. The rapidly rotating 150 star VFTS 682 (12) and the single 90 star 
30Dor 016 (34), ware recently demonstrated to once have belonged to the cluster. The 
proper motion and projected distance from R136 indicates that VFTS 682 and 30Dor 016 
were ejected 29[pc]/30[km/s] 1.0 Myr and 120[pc]/85[km/s] 1.4 Myr ago, respec- 
tively, shortly after the birth of the ^ 2 Myr old cluster. By this time R136 was much too 
young to have experienced a supernova, challenging the favorite explanation for the origin 
of OB runaway stars by the supernova explosion of a stellar companion (3). 

The binary RMS is sufficiently hard to be responsible for having ejected VFTS 682 and 
30Dor 016 from R136, but so massive (M5 ^ 400 Mq (21, 35)) that conservation of linear 
momentum would not impart a high enough velocity kick to its center of mass to eject itself 
from Rl 36. 

Diagonally opposite RMS with respect to the center of R136 is the x-ray source 16, 
which is an ACIS (36) source with a visual as well as a 2MASS counterpart (37). Conser- 
vation of linear momentum may have caused both objects to be ejected in opposite direc- 
tions with respect to R136. Based on the respective projected distances of these two objects 
from the star cluster we predict that RMS has ejected itself and the x-ray source 16 from 
R136. The latter then exhibits a single ^ 130 Mq star or a tight binary. 

RMS and the other massive binary in Fig.O WR2S (HD93162) in the Carina region 
(38), are representative of the population formed dynamically during the core collapse of 
the nearby parent cluster, and flung out by a dynamical interaction. Neither are currently 
in the central portion of the nearby clusters, R136 and Trumpler 16, respectively. Both 
clusters must therefore be in a later stage of core collapse where the binary burning has 
stopped at the cost of the ejection of a rich population of runaways, whereas the core of 
Trumpler 14 collapsed only recently. 

Because each cluster generates one runaway-producing binary during core collapse the 
relative fraction of runaways /run is inversely proportional to the mass of the cluster (Fig. [3]). 
A more complete census of OB runaways can be obtained in the disk of the Milky Way in 
which more than 634 runaway stars have been identified (2). The mass function of these 
runaways is consistent with that of our simulations of young core-collapsed star clusters of 
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6300 Mo (Fig.[3]). 
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Figure 3: The relative fraction of runaways /run as a function of stellar mass m. /run 
represents the number of runaway stars with mass m divided by the total number of stars 
in the cluster with the same mass. The squares, up-pointed triangle and bullets give the 
fraction of runaways found to be associated with the star clusters R136, Westerlund 2 and 
the Galactic field population, respectively. The normalization for the clusters is realized by 
counting the number of cluster members of that mass. The line styles are as in Fig.[Tl 

The majority of the galactic OB runaways seem to originate from star clusters that 
experienced core collapse within the first 1 Myr of their existence. The relative fraction of 
field runaways (Fig. [3]), is somewhat higher than in R136. For stars with a mass ^ 8 Mq 
this fraction is consistent with the stellar ejecta produced in relatively low mass ^ 6300 M0 
star clusters that experience core collapse before they turn 1 Myr old, whereas the R136 
results matches the simulations with a mass of 5 x 10^ Mq. Such star clusters can experience 
core collapse within about a million years if they are born highly concentrated. 



6 



Supporting Online Material 



SOM Text 
Figs. SI toS4 
Table SI 

References (39 — 57) 

In the main paper we discuss the ejection of stars through a gravitational interaction 
from a binary system. Here we present the more technical part of our discussion. In par- 
ticular the discussion on three-body interactions in the cluster core and the macroscopic 
dynamics of the star cluster. This combination of microscopic gravitational dynamics to- 
gether with the global evolution of the star cluster allows us to quantify the processes that 
lead to the ejection of massive runaway stars and at the same time qualify these processes 
by means of global structure evolution calculations of the stellar clusters of interest. Our 
focus is generally on nearby young and massive clusters. From an observational point of 
view these clusters are well studied, and form a pivotal role in the understanding of the 
formation of star clusters. We aim our analysis towards populations of young and massive 
clusters from Tab. 2 of (74). 

The cluster R136 is of particular interest to us, because it is well studied and a rich 
population of runaways was recently identified to roam the neighborhood. With a mass of 
M 60, 000 Mq, a core radius of rcore ^ O.lOpc (13, 14) and a virial radius Tyir — 2.89 pc 
{39) the density profile of R136 fits a King (40) model with a dimension-less central po- 
tential 9.9. As a consequence the density in the cluster core pcore ^ 47, 000 M0/pc^, 
which is quite typical for a core collapsed star cluster (41) (see also Fig. [5]). Inside this clus- 
ter close encounters between stars occur quite frequently. Most common are the encounters 
between a single star and a binary. In §|Al we discuss the results of a series of numerical 
scattering experiments between a binary and a single star. We attempted to match the ini- 
tial conditions for the scattering experiments to mimic the environment in the core of R136 
and in particular the runaway 30Dor 016. In §[B]we report on N-body simulations of star 
clusters which bracket the range of initial masses of observed star clusters, and the cluster 
size is chosen such that they experience core collapse within about one million years. 

A Numerical three-body scattering experiments 

Here we specifically address the relation between the interaction cross section and the 
binary parameters in order to understand the orbital characteristics of the binaries that 
are likely to eject massive stars from a star cluster. The conclusions drawn from the re- 
sults presented in this supplement support the claim in the main paper that relatively wide 
(P 1000 days) and massive binaries are most efficient in ejecting stars from the cluster, 
whereas tight ( ^ 10 day) binaries are unlikely to interact and therefore likely to be primor- 
dial. If binaries with orbital parameters similar to the observed population of short period 
binaries in young and massive star cluster interact dynamically they either cause a collision 
between stars or eject themselves from the cluster. 

We determine the interaction cross section by performing a series of three-body scatter- 
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ing simulations using the sigmaS package within the starlab software environment (42). 
In s igmaS cross sections are determined by executing a number of three-body simulations 
in which a single star is injected into a specific binary. The encounter is initialized using the 
masses and radii of all the three stars, and the relative velocity between the binary center of 
mass and the intruder. The direction from which the intruding star approaches the binary 
and the binary eccentricity are selected randomly. The latter from the thermal distribution 
ranging from a circular orbit to a maximum eccentricity such that the stars do not touch 
at pericenter. The radii of the stars are calculated adopting zero-age main-sequence stars 
with solar composition using MESA (43) called from the AMUSE framework (44). The 
orbital phase and the direction from which the intruder approaches the binary are properly 
randomized (45). The three-body cross section is subsequently determined, as described 
in (46). All calculations are run with a relative energy error AE/E < 10"^. 
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Figure 4: The cross section for an encounter between a single star and a binary as a function 
of the orbital period. The left panel gives the results for single star and binary stars of 
16 M0, whereas the right-hand panel is for 90 stars. The top (blue) curve in both 
panels indicates the cross section for an exchange interaction or a collision between two 
of the stars or all three stars. The green curve gives the cross section for collisions, which 
is a subset of the top (blue) curve. The red curve with bullets gives the cross section for 
encounters that lead to the ejection of one of the stars to a velocity > 30 km/s. The bottom 
two curves give the cross section for encounters that lead to the ejection of a merger product 
(magenta) or of a binary (light blue). Along the right-hand vertical axis the cross section is 
expressed in terms of the encounter rate F for which we adopted cluster parameters similar 
to R136. Along the top horizontal axis we express the encountering binary in terms of 
binding energy in units of kT adopting a mean mass (m) =3.1 M©, which is consistent 
with the average stellar mass in our N-body simulations A, C and D (see §[B]and Tab.[I]). 

In Fig.m we present the results of a series of scattering experiments in which we launch 
a main-sequence star into a binary consisting of two stars of the same mass as the incoming 
star. The left panel in Fig.|4]is calculated using 16 Mq stars with a radius of 4.9 whereas 
the right panel is calculated using 90 Mq stars with a radius of 12.8 Rq. We deliberately 
selected the single stars and the binary components to have the same mass to make the 
interpretation easier, and to prevent selection effects due to the sampling of the other binary 
parameters. We performed additional simulations using stellar masses of 150 200 Mq 
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and 400 The blue dashed curve in Fig.l of the main paper is a least squares fit to the 
stellar mass and the maximum cross section for ejecting a star with a velocity of > 30 km/s, 
which we calculated by means of scattering experiments. This curve is consistent with a 
constant cross section through mass and orbital period. 

The bottom x-axis and left-hand y-axis in Fig. |4] give initial orbital period of the target 
binary (in days) and the cross section (in units of AU^) resulting from the scattering experi- 
ments. Each bullet point in Fig. |4] results from a encounter density of 3000 per surface area 
of the binary, totaling to about 10^ scattering experiments. The various curves are linear 
interpolations between two bracketing points. The top (blue) curve in Fig. |4] gives the total 
cross section for exchanges and collisions (between two or all three stars). The green curve 
presents the cross section for collisions only. After the encounter is resolved, we record 
the velocity of the escaping star. The red bullets give the cross section for an encounter 
which leads to an escaper with an ejection velocity of at least 30 km/s, which is sufficient 
to escape from a cluster like R136 and to be considered a runaway star. The blue ("-h" 
signs) and magenta ("x" signs) curves present the cross section for a merger product and a 
binary with an ejection velocity > 30 km/s, respectively. 

Along the top x-axis of Fig.lH we present the hardness of the binary (in terms of kT, 
where kT is the thermodynamic unit of kinetic energy in the cluster with ^NkT being the 
total binding energy for a cluster of stars) with an orbital separation consistent with the 
bottom X-axis. The binding energy of the binary is calculated adopting cluster parameters 
consistent for R136. Along the right-hand y-axis we express the cross section a (left hand 
y-axis) in terms of the encounter rate between the binary and any other star F (in units of 
Myr~^), which we calculated by adopting a central velocity dispersion of ^disp = 10 km/s. 

A binary of two 16 stars (the minimum for a main-sequence spectral type O star) 
with a ^ 31 year orbital period can already eject an incoming 16 star with a space 
velocity of ^ej — 30 km/s. In the core of a typical young and dense star cluster a binary 
with these parameters exhibits a binding energy of several lOkT. Although such a binary 
is wide, the cross section for producing a runaway in an interaction is relatively small, but 
rapidly increases with decreasing orbital period, until a maximum is reached at a period of a 
few years (see Fig.|4l). The maximum cross section for this binary for producing a runaway 
with ^ej > 30 km/s peaks at an orbital period of 800 days (~ 300 kT) with an encounter 
rate of F 0.01 Myr~^. For higher mass stars (90 Mq) the peak shifts to 2100 days for an 
encounter rate of F :^ 0.3 Myr"^. These rates are consistent with our global estimate in the 
main paper for the orbital separation and mass of a binary that most efficiently ejects an 
incoming star from the cluster. 

For larger orbital period (softer binaries) the cross section for producing a runaway 
drops. For short orbital period (P ^ 100 day) the cross section is dominated by collisions, 
which is consistent with the result presented by (29), although they adopted a = 55 
with stellar masses comparable to those adopted in the bottom panel of Fig.|4l 

The cross section for a collision exceeds the cross section for producing a runaway 
over the entire range of orbital separations, and this difference becomes more prominent 
for wider binaries. The majority of these collision products, remain in the cluster, but a 
small fraction of about 1/10 escape. The fraction of merger products among the escaping 
stars hardly depends on the orbital separation, at least upio P ^ 100 day (compare the red 
buUeted curve with the lower magenta curve in Fig.®. The fraction of collision products 
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among the runaways is then expected to be rather constant, in particular because this trend 
appears to be rather insensitive to the adopted initial conditions. For each five escaping 
normal stars we then expect about double that number of mergers to remain in the cluster 
and one merger product to escape as a runaway. 

The cross section for escaping wide binaries is comparable to that for escaping merger 
products. This cross section does not include < 10 day primordial binaries which dynami- 
cally behave similar fashion as single stars (10). 

B Simulations of massive young star clusters 

To quantify the arguments concerning the orbital characteristics of the binaries that form 
during core collapse, their efficiency of ejecting stars with a high velocity and the distribu- 
tion of stellar masses ejected, we performed a series of direct N-body simulations of star 
clusters over a range of cluster mass and size. 

B.l Method and Initial Conditions 



Name 




N 

^ 'run 


M 


(m) 


Wo 


Tyir 


a 


^rh 


^runaways 








[IO^Mq] 


[Mq] 




[pc] 


[km/s] 


[Myr] 


all m>8 Mq 


A 


2048 


7 


6.4 


3.1 


2 


0.11 


11 


0.37 


18.1 ±2.4 5.0 ±1.3 


B 


8192 


4 


5.7 


0.69 


2 


0.11 


11 


1.2 


42 ±13 3.8 ±1.9 


C 


16384 


6 


51 


3.1 


6 


0.40 


17 


4.4 


27.7 ±5.6 6.2 ±1.0 


D 


65536 


2 


200 


3.1 


8 


1.2 


19 


44 


24.4 ±2.0 3.0 ±1.0 



Table 1: Simulated cluster models. From left to right we identify the model name, number 
of stars in the simulations, the number of runs, the total mass of the cluster, the mean mass 
of the stars, the dimension less depth of the central potential and the virial radius. The 
subsequent columns give derived parameters, which are the central velocity dispersion and 
the half-mass relaxation time. The last two columns give the number of runaways (averaged 
over the A^run performed runs) for all stars which acquire a velocity ^ej > 30 km/s, and for 
the stars with m > 8 Mq that escape. 

In our A^-body simulations we varied the total cluster mass, its size and the stellar mass 
function in order to study the characteristics of the population of runaway stars. In Tab. [J] 
we present an overview of the simulations performed. 

For all simulations we adopted the Salpeter initial mass function (47) between a lower 
limit mmin, which we varied between 0.2 Mq and 1 Mq and a fixed upper limit of 100 Mq. 
The upper limit and the slope of the selected power-law mass function resulted in most of 
the low-mass (M < 10000 M©) clusters to have a maximum mass that rarely exceeded 
30 M0. 

In our simulated star clusters we varied the mass, the virial radius and the concentra- 
tion of the density profile, in such a way that the central density remained roughly the 
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same. Each cluster was initialized using a King (40) model density distribution with a 
dimension-less central potential, Wq, of 2 for the low-mass clusters, 6 for the intermediate 
mass clusters and 8 for the most massive clusters. As a consequence, the central density 
for all clusters - 2 x IO'^MqPc-^ 

In the table (Tab. [I]) we present the initial conditions and main results of our model 
simulations. 
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Figure 5: The evolution of the half-mass density for simulations A, B and C (solid curve). 
Over-plotted using open stars symbols (to the left) are the initial density for the globular 
clusters NGC6397 (49), M4 (48), and 47Tuc (50), The bullets with error bars are densities 
within the effective radius taken from (57). The half-mass density for R136 was taken 
from (13), 

In Fig. [5] we present the evolution of the density within the half-mass radius for one 
of our simulations (C, see Tab. [I]) and for comparison present the initial densities for three 
globular clusters, and the observed densities of several nearby young and massive clusters. 
The half-mass densities for the three globular clusters have been derived by iteratively 
performed A^-body simulations with primordial binaries and stellar evolution until the final 
simulation results matches the observed star clusters (48-50), The initial density used in our 
model is bracketed by those derived for the three globular clusters, and is quite consistent 
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with the observed values. 

We run the simulations using a sixth-order Hermite predictor-corrector integrator (57), 
with a dimension-less accuracy parameter ofr] = 0.2-0.3 (57). We incorporated stellar 
collisions using the sticky-sphere approach for which we adopted the stars to have a radius 
from solar composition zero-age main-sequence stars (52). Initially our simulations have 
no stars with a mass > 100 M©, but due to collisions more massive stars started to appear 
shortly after core collapse. For stars of m > 100 M©, we adopted a mass loss by stellar 
wind of 5.0 x 10~^m [yr~^] (53), for lower mass stars we ignored the mass loss. 

The simulations were curried out using Cray XT4 in National Astronomical Observa- 
tory of Japan with 64-512 cores for 16k and 64k clusters and Little Green Machine in 
Leiden Observatory with 8 cores for 2k clusters. 

B.2 Results from the A^-body simulations 

We ran each simulation up to an age of 3 Myr which is, for all the initial conditions pre- 
sented in Tab.[Il well beyond the moment of core collapse. Core collapse in our models is 
a necessity to develop a mass-segregated mass function and the dynamically active bully 
binary in the core. For each simulation we record the moment of core collapse and the time 
of ejection, mass, and velocity of the escaping stars. 

In Fig. 1 (of the main paper) we present the mass function of the stars in the cluster core 
at an age of 3 Myr. Except for the natural run-to-run variations, the range of initial condi- 
tions did not show a significant difference in the runaway or the core mass function. We 
demonstrate this by presenting the core mass function for all runs combined, but present the 
separate runaway mass function for the three models A, C and D (see Tab.[T]). We compare 
this mass function with the distribution of stellar masses that escaped the cluster, with a 
Selpeter mass function and with the recently observed distribution of massive runaway and 
field stars in the small Magellanic cloud (16). 

The distribution of relatively low mass runaway stars ^ 8 Mq is consistent with the 
Salpeter slope. More massive stars (m ^ 8 Mq) are significantly over represented among 
the runaways and among the population in the cluster core. The mass function of the 
runaway stars with a mass m ^ 8 Mq is indistinguishable from the distribution of masses 
in the cluster center. Although the massive stars are over-represented among the runaways 
and field stars (see Fig. 1 of the main paper), the high-mass end (m ^ 30 Mq) distribution 
is considerably steeper, as was also observed in the small Magellanic cloud (16). The 
consistency between the observed and simulation results further support our finding that 
the majority of massive stars are formed in relatively small (in size as well as in mass) star 
clusters and that the ejection mechanism is dynamic. 

As we discussed in the main paper, both distributions are statistical indistinguishable 
because the cross-section for gravitational interactions is dominated by the relatively wide 
binary formed dynamically during the collapse of the cluster core; the distribution of 
ejected stars depends only weakly on the mass of the single stars. The mass function 
of the runaway stars depends slightly on the mass of the cluster in the sense that the least 
massive clusters tend to have a slightly steeper mass function. This is caused by the lower 
probability of each individual low-mass cluster to contains a massive star. 
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B.3 Velocity distribution 




Figure 6: The cumulative velocity distribution of runaway stars for our simulations with a 
minimum mass of 1 Mq and a total cluster mass of 6300 Mq (left, model A see also Tab.[I]) 
and 5.1 x 10^ Mq (right, model C). The top (solid red) curve gives the distribution of the 
velocities of all runaways with ^gj > 20km/s. For the more massive clusters, the velocity 
distribution of the stars m > 8 Mq (blue dotted curve) and for stars m > 30 Mq (green 
dashes) are somewhat shallower. 



In figure [6] we present the velocity distribution of runaway stars of our simulations with 
a minimum mass of 1 using 6300 Mq and 5.1 x 10^ M©. Opting for 0.2 as a 
minimum mass to the initial mass function does not affect the results significantly. 

The velocity distribution of ejected stars is skewed towards the cluster escape speed, but 
when the binary becomes harder and the ejected star less massive, higher escape velocities 
become more common. The number of runaway star depends on the number of binaries 
formed during core collapse, and not on the mass of the cluster. The distribution of stars 
ejected reflects the hardening process of a single tailor-made binary that formed during core 
collapse. Because it is a single binary that produces all runaways, the scattering process 
in the cluster core drives the energy and rate of escapers. The mass distribution of the 
escapers is then consistent with the stellar population in the cluster core. In a single cluster, 
only one hard binary forms, and the number of ejected stars reflects its hardening process. 
The number of runaways produced, their distribution in mass and in velocity then does not 
depend on the mass of the cluster. 

The characteristic runaway velocity of massive stars, however, does depend on cluster 
mass, even though it is still the single binary process that drives the runaways. The more 
massive stars > 8 Mq tend to have a higher velocity in massive clusters than in the lower 
mass clusters. This trend continues also for a more massive selection of stars, of > 30 M©. 
The lower velocities, on average, of the more massive stars in the lower-mass clusters is 
caused by the smaller proportion of massive stars. Ejecting a massive star from a cluster 
is achieved most easily by a binary consisting of massive stars. However, the lower mass 
clusters have a smaller probability of containing a sufficiently large population of massive 
stars to achieve this. The more massive cluster, therefore, are able to eject massive stars 
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with higher speed. A similar effect was noted in relation to stellar collisions (54). 




Figure 7: The mean escape velocity as a function of mass. We recorded the masses and 
velocities of the stars that have escaped our simulated cluster (models C and D, see Sup- 
plement § B) at an age of 3 Myr. The blue curve is constructed by sorting all escaping stars 
in mass from least to most massive, and measure the means mass and their velocity in a co- 
moving window of 30 stars. We plotted in green the size of the mass-window at the bottom 
of the panel for m = 5 M©, 10 Mq and 20 Mq (squares). The observed mean velocities of 
single main-sequence stars with spectral type A, B and O in the SMC (55) are presented 
in red (top horizontal bars), and the Galactic population of runaways from the Hipparcos 
database (56) is presented in black at 24.4 km/s. 

In Fig.|7] we present the space velocity as a function of stellar mass. The dark blue 
solid curve gives the mean runaway velocity for the simulations C and D as a function of 
mass in a co-moving bin of 30 stars. The horizontal bars (red) give the observed mean 
velocity of stars in the small Magellanic cloud (55). For spectral type O stars (m ^ 16 Mq) 
our simulations are quite consistent with the observed velocities, although for later type 
spectral types our simulations tends to somewhat under produce the velocities. Although 
using the Hipparcos catalogue (56) have derived a mean velocity of 24.4 ± 1.2 km/s for 
runaways with a mean mass of 5.4 ± 4.5 Mq. The average runaway velocity for stars with 
m ^ 16 Mq seem to be systematically higher in the SMC than those observed in the Galaxy 
and our simulations. 
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